Berezinskii-Kosterlitz-Thouless transition in two-dimensional dipole systems 
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The superfluid to normal fluid transition of dipolar bosons in two dimensions is studied throughout 
the whole density range using path integral Monte Carlo simulations and summarized in the phase 
diagram. While at low densities, we find good agreement with the universal results depending only 
on the scattering length a s , at moderate and high densities, the transition temperature is strongly 
affected by interactions and the elementary excitation spectrum. The results are expected to be of 
relevance to dipolar atomic and molecular systems and indirect excitons in quantum wells. 



Dipolar bosonic systems are of increasing interest for 
various recent experiments studying the onset of su- 
perfluidity in nonideal Bose systems and its connection 
with correlation and quantum degeneracy effects. Exam- 
ples include dipolar gases, as in recent studies of 52 Cr 
atoms p], bosonic molecules, e.g. SrO, RbCs, LiCs, 
40 K 8 'Rb [2], as well as indirect excitons in semiconductor 
quantum wells [3J 0] • A number of theoretical and com- 
putational studies have addressed the properties of two- 
dimensional (2D) repulsive dipolar bosons at zero and 
low temperatures j5rfZ|. They include the ground-state 
energy, the structural and coherence properties, such as 
the one-body density matrix and the condensate fraction. 
Quantum Monte Carlo studies at T = have covered the 
whole range of coupling strengths up to the crystalliza- 
tion transition. However, the finite-temperature proper- 
ties, in particular, the superfluid transition temperature 
T c , remain unexplored. 

Previous numerical investigations for 2D homogeneous 
Bose gases [HI E] have shown that in the dilute (weak 
coupling) regime, na 2 < 1CP 2 , where n is the density 
and a s the s-wave scattering length, the exact shape of 
the interaction potential is irrelevant for T c which is a 
function of no 2 only [TU] 
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where the numerical coefficient is v — 380(3) How- 
ever, for moderate and high densities where correlation 
effects are important no analytical expression is avail- 
able for T c , calling for investigations by direct numerical 
simulations which is the main goal of the present paper. 
By performing first principle path integral Monte Carlo 
(PIMC) simulations we demonstrate that, with increas- 
ing interaction strength, the superfluid phase is first sta- 
bilized (T c increases) and then destabilized and vanishes 
when the system forms a dipolar solid. We present an ex- 
planation of this non-monotonic behavior of T c suggest- 
ing that it arises from a competition between two types 
of elementary excitations in the normal fluid component 
- phonons and rotons. 



Model and parameters. We focus on a pure 
dipole model relevant e.g. for various bosonic atoms or 
molecules and indirect excitons at low densities where 
the dipole moment is a free parameter which can be ex- 
ternally controlled, e.g. via an electric field (UH1]. The 
2D dipole system is described by the Hamiltonian 
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which is brought to a dimensionless form using the units 
a = l/\/n and Eq = h 2 /ma 2 . The system properties are 
defined by the dipole coupling D = p 2 /ei ) a 3 Eo and the 
temperature, T = ksT/EQ. The thermodynamic equi- 
librium states of system Q were sampled by PIMC sim- 
ulations with the worm algorithm |13| . We studied the 
system Q from weak to strong coupling, D = 0.01 ... 20. 
In agreement with Refs. [5HZ] we observe formation of a 
crystalline phase at D ~ 18. 

Superfluid transition and phase diagram. In 2D 
the superfluid-normal phase transition occurs at a finite 
temperature T c and follows the Berezinskii-Kosterlitz- 
Thouless scenario (BKT) induced by interaction ef- 
fects [13] • To obtain a reliable result for T c in the thermo- 
dynamic limit, T c (oo), from simulations of a finite system 
of size L = y/N, we apply a finite-size scaling analysis to 
T C {L). We assume the essential singularity [T^ of the 
correlation length £(T) ~ e at ~ 1/2 ,t = (T/T c - l)\ T -yT e , 
with a being a non-universal temperature-density depen- 
dent scaling factor. Near T c the role of £ is taken over by 
L, leading to (b is a constant, cf. Fig. [T]b) 



T C (L) = T c (oo) 
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where T C (L) is determined by the scenario of the uni- 
versal jump of the superfluid fraction (second equa- 
tion) [16] . Here, the superfluid density n s is obtained 
via the winding number estimator |17| . n s (L,T) = 
mk B T(W 2 )(T,L)/2h 2 which is directly evaluated by 
PIMC simulations. Combining this with Eq. T C (L) 
is determined by the condition (W 2 )(T C (L), L) = 4/ir. 
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FIG. 1: (a) Winding number (W 2 ){T) for 'different system 
sizes: N = 36, 49, 81, 324. Dipole coupling D — 5. Horizontal 
dashed line: W 2 (T C ) = 4/tt. (b) System size dependence of 
the critical temperature for several coupling strengths D. 



In Fig. [T] we show the temperature dependence and the 
finite-size scaling of (W 2 )(T, L) and T C (L). We observe 
that T C (L) shifts systematically with L to lower values 
(Fig.Jlji). The extrapolation to the thermodynamic limit, 
T c (oo), fitting the simulation data by Eq. |3|, is reported 
in Fig. [TJd. For strong coupling (D > 10) the finite-size 
corrections to T c in Eq. ^ become important, therefore, 
we excluded from the fit the smallest system (N = 36). 

Using the extrapolated data, T c (oo,D), we construct 
the phase diagram, cf. Fig. [2] and Table |T] At small cou- 
pling, D < 0.01, our results are well reproduced by the 
asymptotic expression ([TJ, i.e. here, details of the in- 
teraction potential are not important. However, the va- 
lidity range of Eq. ([lj is limited to very low densities, 
na 2 < 10~ 3 [computing the scattering length from the 
solution of the Schrodinger equation for 2D dipoles gives 
the relation na 2 ss 10.05-D 2 ]. This density is an order 
of magnitude smaller than for a 2D gas of hard disks [9_ 
indicating that the long-range character of the dipole in- 
teraction causes substantially earlier deviations from the 
dilute gas limit. 

Now we analyze, in more detail, the change of T c with 
coupling, Fig. [2} For small D, T c monotonically increases 
and reaches a maximum around 1 < D < 2 whereas, for 
larger couplings, T c monotonically decreases again until 
the system freezes into a non-superfluid dipole crystal. 
[This is preceded by a narrow hysteresis region (15 < D < 
18), shown by the vertical dashed lines, where a bubble- 
type structure is expected |19| which, however, is beyond 
the scope of the present paper.] In the following, we 
demonstrate that the non-monotonic behavior of T C (D) 
can be explained by the excitation spectrum. 

Excitation spectrum. The excitation spectrum u)(q) 
can be approximated from above through the finite tem- 
perature generalization of the Feynman relation [20 
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where S(q, T) is the static structure factor which we com- 
puted in the PIMC simulations for a broad range of tem- 
peratures, 0.5 <T< 3.3, below and above T c , cf. Fig. [3^. 
As D approaches the crystallization point, a sharp peak 
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FIG. 2: (a): Phase diagram in the T-D plane. The system 
is superfluid below the solid (blue) line. The two vertical 
dashed lines bound the gas-solid hysteresis region to the right 
of which the system is in a dipole crystal phase 5-7 18J. The 
horizontal dashed line is the upper bound for T c obtained 
by replacing n s — > n. The dotted line denotes the estimate 
for T c according to Eq. (TTT) . (b): D-dependence of the normal 
density n n (PIMC result )7 the quasiparticle contribution n^ p , 
Eq. Q and the phonon and roton contributions, Eqs. Q and 
»|), respectively, for T = T c . The difference between n n and 
equals and is due to vortices, cf. table [i] 



develops near the wave number q corresponding to the 
mean interparticle distance, q a = 2ir. Note that, while 
S(q, T) shows some T-dependence for qa < 3, the Feyn- 
man spectrum, LOF(q,T), stays almost unchanged in a 
broad temperature interval, T < 3.3, and is close to 
the ground state result [24]. Therefore, the spectra 
shown in Fig. [4] for T = 0.5 are representative for the 
low-temperature behavior. 

In the long wavelength limit, qa — > 0, ujp becomes 
exact yielding a linear dispersion, u)p{q) — c s q, cf. Fig. [4^ 
which is in agreement with the result for classical 2D 
dipoles Our results for the sound speed, c s (T) = 

ujf(q> T)/q\ q -yQ, extracted from the data for N = 324 
particles are summarized in Tab. [I] and agree within 4% 
with the ground state values of Ref. [7] . 
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FIG. 3: (a) Static structure factor and (b) density response 
function at T = 0.5 (solid lines) and T = 3.3 (dashed lines) 
for three couplings, D = 0.1, 1, 10 (numbers in the figure). 

A significant improvement of the spectrum is achieved 
by using a sum-rule approach [25, 26 by combining the 
PIMC results for S(q, T) and the static density response 
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x(q,T), yielding a rigorous upper bound 

hoj(q) < hu: x (q, T) = 2nS(q, T)/ X (q, T), (5) 

X(q) = -J F(q,r)dT, F(q,r) = - (p q (r)p q (0)) , 

where x(g) is obtained from the imaginary-time density- 
density correlation function F(q, r) which has been di- 
rectly evaluated in PIMC simulations. With the increase 
of D (Fig. [3]d) the response function sharpens and shifts 
continuously towards q . 

In Fig. [4] we show lu x , Eq. ^ , together with the Feyn- 
man spectrum and the correlated basis functions (CBF) 
result |57| at four dipole couplings. All three approxima- 
tions show the same general trend which resembles su- 
perfluid helium: with increasing coupling the spectrum 
develops a roton minimum at finite q sa q which becomes 
deeper with increasing D. While for qa < 1.5 (sound 
range) all approaches are in quantitative agreement, for 
qa > 2 the Feynman approximation becomes inaccurate. 
Its error increases with D and exceeds 100% at the crys- 
tallization point for uj at the minimum. The PIMC result 
to x (q) agrees suprisingly well with wcbf(<z)- Our simula- 
tions predict a deeper minimum Lu(qa) and are expected 
to be more accurate here. Further, for q > 7.5, the upper 
bound u> x (q) approaches a free-particle spectrum (simi- 
lar to LJp), except for D > 15, whereas CBF, at strong 
coupling (D = 15), shows the onset of a plateau. In anal- 
ogy with superfluid 4 He a plateau might be expected at 
twice the roton minimum energy but it appears that all 
schemes violate this threshold and this question requires 
additional research. 




FIG. 4: Excitation spectrum e q = hui{q, T) for T = 0.5 and 
different couplings D. PIMC results, wp, Eqs. Q, and oj x , 
Eq. uSp , are compared with the CBF spectra (D = 1, 4, 8, 16), 
Ref. |27| and (hq) 2 /2rn (short dashed lines). 



We now return to the phase diagram. T C (D) in Eq. (|3j, 
is related to the normal density by n s {T) = n — n n (T) 



which for non-interacting quasiparticles is determined by 
the excitation spectrum s p — hw x (p; D, T) via the Lan- 
dau formula j28] 
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where n B = [ef^-l]- 1 is the Bose distribution function, 
and we used m = 1. The superscript QP denotes the 
dilute quasiparticle approximation. 

At low temperatures the main contribution to n^ p , 
Eq. ([6]), comes from the lowest energy quasiparticles - 
phonons - with the dispersion e 9 ^ 1 = cp, and rotons, e r v = 



A + (p — po) /2/i, with the approximate result [28 in 2D 
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The free parameters in the dispersions e^ 1 and e r p are 
determined from a fit to the simulation result huj x (p). 
The values for c (c ~ c s ) and the roton gap A are 
summarized in Table |T] We note that A in the liq- 
uid phase decreases exponentially with the dipole cou- 
pling: A(D)/Eq — a\ exp(— a 2 D — a 3 D 2 ), with the 
best fit parameters: ai = 15.11(5), a 2 = 0.088(2), 
a 3 = —0.00120(8), whereas at the crystallization point 
we find A(D = 18) /Eq = 4.57. 

Consider now the D-dependence of nP h and rf n , cf. 
lower part of Fig. [2j At weak coupling, D < 1, there is 
no roton minimum in the spectrum, cf. Fig. [4] D = 1, 
and rf n = 0. At the same time there is a strong phonon 
contribution, nf/ 1 , which monotonically decreases with D 
since, at a fixed temperature, the sound speed increases 
with coupling, c s ~ D 1 / 2 , cf. Tab. [i] As soon as the 
spectrum contains a roton minimum, i.e. for D > 1.5, cf. 
Fig. [4] phonon excitations are practically negligible and 
the rotons dominate which is due to the larger density 
of states (~ qdq). The roton density n r n monotonically 
increases with D, cf. Fig. [2jb due to the reduction of 
the roton gap, cf. Tab. [T] These competing trends of 
phonon and roton excitations give rise to a minimum in 
the sum n,P h (D) + n r n {D) at a coupling around D = 1.5. 
Note that this sum of the approximate contributions ([7| , 
(l8| rather well reproduces the full quasiparticle density 
r% F (D), Eq. §. 

Most interestingly, the position of the minimum is very 
close to the maximum of the superfluid transition tem- 
perature T c which is observed in the range D ~ 1 — 1.75, 
and the T C (D) appears to follow just the opposite be- 
havior of the curve n^ p (Z3). In fact, this opposite trend 
could be expected from Eq. Q stating that T c ~ n — n n 
where n n (T c ) is the full normal density at the critical 
point. However, there is no obvious reason why n n {D) 
should follow n,^ p in the case of strongly interacting 
bosons. Our simulations allow us to directly compare 
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the exact value, n n to n% p . As can be seen in Fig. |2|b 
there is a dramatic difference between the two indicating 
the existence of an additional contribution to the nor- 
mal density, i.e. n n — +«^, which is by far dominant 
in the region around the maximum of T c exceeding 
by as much as a factor thirty. 

The physical origin of are interaction effects and 
vortices missing in the model (J6j> . Large scale vor- 
tices are a key in the BKT theory and they are known 
to substantially reduce the true macroscopic superfluid 
density n s compared to its local value n l ° c (T). It has 
been shown that this suppression is similar to screening 
[TH, |2"9"] characterized by an effective dielectric function, 
n s = n[ oc /e. For the present system, it is reasonable to 
identify n l ° c — > n — n^ p which relates e to rC n according 
to n w n = (n — n^ p )(l — e _1 ). Our simulations allow us 
to directly compute the vortex density and e(D) for all 
couplings, the results are included in table [I] Notice that 
and e rapidly increase for strong coupling. However, 
for D < 3, i.e. in the range where T C (D) reaches its 
maximum they are almost D-independent. This behav- 
ior is constrast to the one of n^'(D) and n r n {D). Thus 
we may conclude that - despite the small absolute value 
of n I ^ l {D) and n r n (D) compared to n v n , it is the phonon 
and roton excitations which are responsible for the shape 
of the phase boundary T C (D) and for the stabilization of 
the superfluid phase around D ~ 1 — 1.75. 

TABLE I: Coupling parameter dependence of the superfluid 
transition temperature T c , the superfluid fraction j 3 (T c ) = 
2mT c /nh 2 n, the sound speed c s [in units of the dipole fre- 
quency, [o;£>a], ~ 27rp 2 n/(ma 3 ), a = (7m) -1 / 2 ], the vortex 
density n^(T c ) and the effective dielectric function, e(T c ). 
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In conclusion, the finite-temperature phase diagram 
of a 2D dipole system has been investigated by first- 
principle PIMC simulations over the entire coupling 
regime. We found that the superfluid density at T c does 
not exceed 90% and drops to about 58% near the crystal- 
lization point. The normal density is dominated by vor- 
tices and contains a small fraction of phonons and rotons. 
Yet the competition of the latter two is responsible for 
the observed non-monotonic behavior of T C (D). Further- 
more, an upper bound for the single-particle spectrum 
has been computed which significantly improves the re- 



sult of the Feynman approximation. We expect that our 
predictions are of direct relevance for experiments with 
atomic and molecular dipole systems as well as for indi- 
rect excitons in semiconductor quantum wells. 
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